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LONG-TERM  GOALS 


£  '  *> 


The  primary  goal  of  this  project  is  to  substantially  reduce  the  computational  complexity  of  the  adap¬ 
tive  beamformer  employed  in  real-time  passive  sonar  systems  for  underwater  surveillance.  The  main 
methodology  for  achieving  this  computational  reduction  is  reduced-rank  adaptive  filtering  via  the 
method  of  Conjugate  Gradients.  Reduced-rank  adaptive  filtering  is  vitally  important  for  underwa¬ 
ter  surveillance  due  to  the  large  number  of  hydrophones  typically  employed,  coupled  with  issues  of 
sample  support  related  to  target  mobility  and  the  high  degree  of  nonstationarity  of  the  underwater 
environment. 

OBJECTIVES 

A  key  objective  is  to  develop  a  computationally  efficient  implementation  of  the  Multi-Stage  Wiener  Fil¬ 
ter  (MWF)  for  adaptive  beamforming  in  scan  mode.  The  Minimum  Variance  (MV)  Spatial  Spectrum 
estimate  for  a  particular  bearing  may  be  computed  as  the  reciprocal  of  the  inner  product  between 
the  “look  direction”  vector  for  that  bearing  angle  and  the  solution  to  the  MVDR  linear  system  of 
equations,  with  the  sample  covariance  matrix  on  the  left-hand-side  and  the  “look  direction”  vector 
on  the  right-hand-side.  Explicitly  computing  the  inverse  of  the  sample  covariance  matrix  is  computa¬ 
tionally  burdensome  and  problematic  if  the  sample  support  is  not  adequate.  The  Multi-Stage  Wiener 
Filter  (MWF)  provides  a  reduced-rank  solution  to  the  Minimum  Variance  Distortionless  Response 
(MVDR)  equations  that  is  computationally  efficient  as  well  as  robust  to  sample  support.  However,  in 
scan  mode,  the  MVDR  equations  need  to  be  solved  for  each  of  a  large  number  of  “look  directions.” 
This  has  limited  the  application  of  MWF  in  MVDR  based  scanning.  However,  the  recently  discovered 
equivalence  between  MWF  and  the  iterative  search  method  of  Conjugate  Gradients  (CG)  provides 
the  key  to  computationally  efficient  scanning  via  the  MWF:  the  solution  to  the  MVDR  equations  for 
one  bearing  angle  is  used  to  initialize  the  CG  search  for  the  solution  to  the  MVDR  equations  for  the 
next  bearing  angle  on  the  grid.  As  a  result,  only  a  few  steps  of  CG  are  needed  each  time  the  “look 
direction”  is  incremented  in  the  process  of  scanning. 

A  second  objective  is  to  develop  efficient  FFT  based  implementations  of  the  CG  method  for  gener¬ 
ating  the  Minimum  Variance  spatial  spectrum  that  avoid  the  need  to  explicitly  estimate  the  covariance 
matrix.  A  third  objective  is  conduct  extensive  simulations  to  verify  the  efficacy  of  this  approach  in 
complex  underwater  environments.  These  simulations  will  also  include  an  alternative  to  CG  based  on 
Constrained  Steepest  Descent.  A  fourth  objective  is  to  use  both  theoretical  sensitivity  analysis  and 
the  results  of  these  simulations  to  develop  rules  for  variable  grid  step  size  and  the  number  of  CG  steps 
at  each  grid  point. 

APPROACH 

The  search  volume  of  interest  is  scanned  by  discretizing  the  bearing,  range,  and  depth  location 
variables.  For  each  “look  direction”  of  interest,  one  ostensibly  has  to  solve  the  MVDR  equations 
Rxiw(0)  =  s (0),  where  Rzz  is  the  sample  covariance  matrix,  and  s(0)  and  w(0)  are  the  steering 
vector  and  adaptive  beamforming  weight  vector  for  the  particular  look  direction  0.  In  scan  mode,  the 
MVDR  equations  need  to  be  solved  for  each  of  a  large  number  of  look  directions.  For  each  prospective 
look  direction,  MWF  ostensibly  starts  “over  again”  by  using  the  current  look  direction  vector  s(0)  as 
the  first  basis  vector  for  the  forward  recursion.  Such  a  mode  of  operation  has  limited  the  application 
of  MWF  in  MVDR  based  scanning. 

Yet,  a  primary  advantage  of  MWF  deals  with  converging  to  the  near-optimal  adaptive  filter  during 
a  period  when  the  underlying  signal  statistics  are  relatively  stationary.  The  larger  the  number  of 
elements  comprising  the  array,  the  longer  it  takes  for  the  adaptive  beamformer  to  converge  to  a  point 
where  it  is  providing  a  near-optimal  SINR,  for  example.  If  the  underlying  signal  statistics  are  varying 
too  quickly,  a  high-dimensional  beamformer  may  not  converge  and  thus  lose  track.  MWF  based 
reduced-rank  adaptive  filtering  offers  substantial  decreases  in  convergence  time,  leading  to  enhanced 
performance  in  nonstationary  target  scenarios. 
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Fortunately,  the  CG  connection  provides  the  key  to  computationally  efficient  scanning  via  the  MWF. 
CG  is  an  iterative  search  scheme;  CG  takes  less  steps  to  converge  to  the  “desired  solution”  the  closer 
the  starting  point  is  to  the  “desired  solution”.  The  solution  to  RxXw(0)  =  s(0)  should  be  used  to 
initialize  the  CG  based  search  for  the  solution  to  Rxxw(0  +  <5)  =  s(0  +  <5),  If  S  is  small  enough,  only 
two  or  three  steps  of  CG  are  adequate  for  each  new  look  direction.  The  first  look  direction  will  have  to 
be  solved  in  the  “standard”  MWF  manner,  with  the  first  basis  vector  of  the  forward  recursion  equal 
to  the  steering  vector  for  the  initial  “look  direction.” 

The  other  key  implication  of  the  CG  discovery  is  that  MWF  may  be  used  in  conjunction  with  Principal 
Components  (PC)  based  rank  reduction.  Being  an  iterative  search  method,  CG  can  be  applied  in  a 
lower-dimensional  space,  the  signal  subspace,  for  example.  The  beamformer  for  each  target  being 
tracked  needs  to  adapt  with  time.  Dr.  Norm  Owlsey  and  others  have  shown  that  the  use  of  Principal 
Components  lowers  the  sample  support  needed  for  estimating  the  covariance  matrix  from  roughly 
twice  the  number  of  sonar  sensors  to  twice  the  number  of  sources.  Coupled  with  computationally 
efficient  algorithms  for  tracking  the  signal  subspace,  PC  can  effect  a  substantial  reduction  in  the 
computational  load  by  effectively  compressing  the  high-dimensional  “element  space”  covariance  matrix 
to  a  covariance  matrix  of  dimension  equal  to  the  number  of  sources,  denoted  Rcc.  CG  may  be  applied 
as  follows:  the  solution  to  Rcc(n)w(0(n))  =  s(9(n))  is  used  to  initialize  the  CG  search  for  the  solution 
to  Rcc(n  4-  D)w(0(n  +  D))  =  s(6(n  +  D)),  where  Rcc(n)  and  Rcc(n  +  D)  are  the  compressed  sample 
covariance  matrices  at  time  n  and  n  +  D,  respectively,  and  9(n)  and  9(n  +  D)  are  the  look  directions 
for  the  target  of  interest  at  time  n  and  n  +  D ,  respectively.  The  value  of  D  depends  on  the  update 
rate. 

WORK  COMPLETED 

Key  accomplishments  during  the  first  year  effort  were  fueled  by  our  recent  discovery  of  an  inherent 
relationship  between  the  MWF  and  the  Iterative  Search  Method  of  Conjugate-Gradients  (CG).  This 
lead  to  a  host  of  new  algorithmic  implementations  of  CG-MWF  that  offer  a  number  of  important 
advantages  over  the  original  implementation  of  MWF.  The  discovery  that  adding  a“stage”  to  the 
MWF  was  equivalent  to  taking  a  “step”  of  a  CG  search  has  lead  to  substantial  improvements  to 
the  MWF  in  terms  of  (i)  computationally  efficient  implementations  that  exploit  structure  in  the 
data  matrices  or  correlation  matrices,  such  as  Toeplitz,  sparseness,  etc.,  (ii)  amenability  to  real-time 
implementation,  (iii)  amenability  to  “smart”  initialization  (based  on  information  learned  or  statistics 
estimated  during  the  process  of  searching  for  a  training  sequence  embedded  in  the  data,  for  example), 
(iv)  easy  incorporation  of  a-priori  information,  constraints,  and/or  “past  history”  for  accelerated 
convergence,  and  (v)  amenability  to  complementary  PC  based  rank  reduction. 

The  basic  implementation  of  CG-MWF  for  solving  the  N  x  N  system  of  Wiener-Hopf  equations 
RxzW  =  Tdx  is  summarized  in  Table  1,  where  r^x  is  the  cross-correlation  vector  between  the  desired 
signal  and  the  observed  data.  For  the  beamforming  application,  r^x  is  the  steering  vector  for  the 
current  look  direction.  We  have  formulated  new  algorithmic  implementations  of  CG-MWF  that  offer 
a  number  of  important  advantages  over  the  original  implementation  of  MWF.  In  addition,  CG-MWF 
offers  a  number  of  VIP  advantages  over  both  LMS  and  RLS  based  adaptive  filtering. 

The  basic  CG-MWF  algorithm  is  delineated  on  the  left  hand  side  of  Table  1.  Each  step/iteration  of 
CG-MWF  involves  two  vector-vector  multiplies,  three  simple  vector  updates,  and  one  matrix-vector 
product.  The  latter  is  the  most  computationally  burdensome  step  of  CG-MWF,  Rxx  may  be  estimated 
from  M  blocks  of  data  as  Rxx  =  XXH  /M,  where  X  is  an  N  x  M  Toeplitz  data  matrix.  Through 
the  circulant  extension  of  a  Toeplitz  matrix,  the  code  in  the  right  hand  side  of  Table  1  shows  that  the 
product  R^Uj  can  be  computed  via  FFT’s  leading  to  two  VIP  advantages:  (1)  substantially  reduced 
computational  complexity  by  choosing  the  FFT  length  to  be  a  power  of  two  and  (2)  substantially 
reduced  memory  requirements  as  Rxx  =  XXH/M  need  not  be  formed  and  thus  not  stored.  Relative 
to  the  latter  point,  all  that  is  required  is  a  one-time  FFT  of  the  data  stream  from  which  the  M 
successive  data  blocks  are  implicitly  extracted. 
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Recent  algorithmic  implementations  of  MWF  by  Ricks  and  Goldstein  also  avoid  having  to  form  (or 
store)  a  covariance  matrix.  In  contrast  to  the  FFT  based  CG-MWF  algorithm  described  above,  the 
innovative  lattice,  modular  implementation  of  MWF  of  Ricks  and  Goldstein  requires  a  backwards 
recursion  as  well  as  a  forward  recursion  and  it  does  not  directly  yield  the  Wiener-Hopf  weights.  The 
relationship  between  the  standard  version  of  MWF  and  their  innovative  lattice/modular  implemen¬ 
tation  of  MWF  is  analogous  to  the  relationship  between  the  standard  version  of  RLS  and  fast  RLS 
algorithms. 

As  part  of  the  second  year  effort,  we  are  developing  further  refinements  to  CG-MWF  that  advan¬ 
tageously  exploit  the  fact  that  it  works  directly  on  the  covariance  matrix.  In  contrast,  RLS  either 
explicitly  or  implicitly  works  with  the  inverse  of  the  covariance  matrix.  As  a  result  of  the  nonlin¬ 
ear  relationship  between  the  covariance  matrix  and  its  inverse,  the  innovations  we  propose  have  no 
counterpart  in  the  realm  of  RLS  based  algorithms.  As  a  result,  in  addition  to  the  huge  advantage 
of  CG-MWF  over  RLS  under  conditions  of  low  sample  support,  due  to  the  inherent  reduced-rank 
processing  of  CG-MWF,  there  are  also  major  advantages  in  terms  of  initialization  and  computational 
complexity.  The  computational  complexity  reduction  of  CG-MWF  relative  to  RLS  arises  primarily 
from  the  former’s  ability  to  exploit  the  Toeplitz  structure  of  the  covariance  matrix  via  FFT’s.  Since 
the  inverse  of  a  Toeplitz  matrix  is  not  Toeplitz  (in  general),  RLS  does  not  offer  this  feature. 

Another  key  advantageous  feature  of  CG-MWF  relative  to  RLS  is  “smart  initialization.”  Equalization 
of  underwater  digital  communication  signals  -  equalizing  distortion  arising  from  propagation  effects 
-  is  used  in  the  illustrative  simulation  example.  In  this  application,  one  may  obtain  an  initial  esti¬ 
mate  of  the  impulse  response  of  the  propagation  channel  through  semi-blind  means  (limited  training 
plus  exploitation  of  signal  properties.)  We  say  “initial”,  since  the  channel  is  time- varying,  in  general. 
The  initial  channel  estimate  may  be  used  to  initialize  the  covariance  matrix,  R^,  in  the  CG-MWF 
algorithm,  as  well  as  the  cross-correlation  vector,  r dx-  The  former  is  the  outer  product  of  the  con¬ 
volution  matrix  of  the  channel  with  itself,  while  the  latter  is  is  the  channel  itself  with  some  zero 
padding.  Both  quantities  are  then  updated  on  a  per  snapshot  basis  as  time  evolves  to  implicitly  track 
the  time-variations  of  the  channel.  Preliminary  results  indicate  that  this  channel  based  initialization 
of  CG-MWF,  which  is  possible  due  to  the  fact  that  CG-MWF  works  directly  with  Ria;,  has  a  very 
dramatic  effect  on  the  convergence  speed  of  CG-MWF.  In  contrast,  even  if  an  estimate  of  the  channel 
is  available,  there  are  no  means  through  which  RLS  can  utilize  it  since  that  would  require  inverting 
the  outer  product  of  the  convolution  matrix  of  the  channel  with  itself;  this  totally  defeats  the  purpose 
of  RLS  in  addition  to  representing  a  huge  attendant  computational  burden. 

Although  the  target  application  is  underwater  surveillance,  the  CG-MWF  algorithm  is  as  universally 
applicable  as  the  classical  adaptive  filtering  schemes  of  LMS  and  RLS.  In  fact,  I  believe  that  CG 
represents  the  next  generation  of  adaptive  filtering  algorithms. 

RESULTS 

A  simple  illustrative  simulation  is  presented  to  demonstrate  the  efficacy  of  the  efficient  FFT  based 
CG-MWF  algorithm.  The  application  is  linear  equalization  of  a  multipath  distorted  QPSK  signal.  In 
addition  to  a  direct  path,  there  is  a  70%  ghost  at  a  half-symbol  delay  with  a  phase  of  165°.  Factoring 
in  the  pulse  shaping,  the  effective  channel  length  was  roughly  L  =  13.  An  equalizer  length  of  Ng  =  20 
was  selected;  the  taps  were  symbol-spaced.  N  =  32  is  the  FFT  length. 

The  multipath  channel  is  frequency  selective,  thereby  causing  the  intersymbol  interference  shown  in  the 
received  constellation  plotted  in  the  upper-left  corner  of  Figure  1.  Applying  the  respective  equalizing 
weight  vectors  obtained  after  3  steps  and  5  steps  of  FFT  based  CG-MWF  yields  the  constellations 
plotted  in  the  lower-left  and  lower-right  corners  of  Figure  1,  respectively.  This  demonstrates  that 
near-optimal  performance  may  be  obtained  with  a  few  steps  of  CG-MWF,  and  that  CG-MWF  may 
be  implemented  efficiently  via  FFT  processing.  In  addition  to  the  savings  associated  with  the  FFT 
based  computation  of  the  elements  of  the  matrix- vector  product  RxxUi  at  each  step  of  CG,  we  also 
avoid  the  computation  associated  with  explicitly  estimating  Rxx. 
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IMPACT/APPLICATIONS 


This  work  will  yield  substantial  reductions  in  the  computational  complexity  associated  with  imple¬ 
menting  adaptive  beamformers  in  real-time  passive  sonar  systems  for  underwater  surveillance. 

TRANSITIONS 

An  all-day  tutorial  entitled  “TUTORIAL  ON  REDUCED-RANK  ADAPTIVE  FILTERING  BASED 
ON  THE  MULTI-STAGE  WIENER  FILTER”  was  presented  on  18  June  2002  to  the  CITE  Group  at 
Rome  Labs  headed  by  Dr.  Bruce  Suter.  Dr.  Bruce  Suter  and  CITE  will  implement  and  assess  the 
performance  of  the  Conjugate-Gradient  based  reduced  rank  filtering  schemes  developed  as  part  of  this 
project.  We  will  continue  to  report  of  this  ongoing  transition  as  it  progresses. 

RELATED  PROJECTS 

There  is  a  synergism  between  this  work  and  research  being  conducted  in  parallel  for  a  National 
Science  Foundation  project  entitled  “Reduced-Dimension  Decision  Feedback  Equalizers  (DFE’s)  for 
4G  High-Speed  Wireless  Digital  Communications.”  This  is  a  single  Principal  Investigator  grant  funded 
by  the  Computing  and  Communications  Research  Division  of  the  CISE  Directorate  at  NSF.  (Grant 
Number:  CCR-0118842.  Duration:  1  Sept.  2001  -  31  Aug.  2004.)  Equalization  of  Digital  TV  is  an 
ideal  application  for  CG  based  adaptive  filtering  due  to  the  high  dimensionality  of  the  DFE;  current 
receiver  chips  have  576  DFE  taps  that  need  to  be  adapted  with  time.  In  addition  to  the  synergism 
relative  to  theoretical  developments,  there  is  a  synergism  emerging  with  respect  to  the  underwater 
application  following  from  a  lecture  recently  delivered  at  the  Second  Sensor  Array  and  Multichannel 
Workshop  held  4-6  August  2002  in  Rosslyn,  VA.  The  lecture  was  by  Dr.  Johann  F.  Bhme  of  the 
Ruhr-University  in  Bochum,  Germany,  and  was  entitled  “Multichannel  Blind  Equalization  in  Digital 
Underwater  Acoustic  Communications:  Experimental  Results.” 

PUBLICATIONS 

•  Michael  Zoltowski,  “Conjugate  Gradient  Based  Adaptive  Filtering  with  Application  to  Space- 
Time  Processing  for  Wireless  Communications,”  Second  IEEE  Sensor  Array  and  Multichannel 
Signal  Processing  Workshop  (SAM  2002),  Rosslyn,  VA,  5-6  August  2002. 

•  Matthew  Weippert,  John  Hiemstra,  J.  Scott  Goldstein,  and  Michael  Zoltowski,  “Insights  from 
the  Relationship  Between  the  Multistage  Wiener  Filter  and  the  Method  of  Conjugate  Gradients,” 
Second  IEEE  Sensor  Array  and  Multichannel  Signal  Processing  Workshop  (SAM  2002),  Rosslyn, 
VA,  5-6  August  2002. 
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L=block  length;  N=FFT  length; 

M= weight  vector  length;  N=M+L-1 
X=fft(xd,N); 

F=(exp(-j*2*pi/N).fO:N-l]’).*conj(X); 
w=zeros(M,l);  u=rdx;  g=-u; 
i=g’*g; 

for  i=l:Nstop, 

d=ifft(X.*fft(u,N),N);  y=d(end-L+l:end!l); 
z=ifft (F.*fft (y,N) ,N) ;  v=z(end-M+l :end,l ) ; 
eta=l/(u’*v); 
w_old=w; 
w=w_old+eta*u; 
g.old=g; 
g=g_old+eta*v; 

Lold=l; 
i=g!*g; 
psi=l/l_old; 
uold=u; 

u=-g+psi*u_old; 
end 

Table  1.  Standard  CG-M\¥F  (left)  and  FFT  Implementation  of  CG-MWF  (right) 


Direct  Block  CG-MWF. 


Recvd  Signal  Constellation  Computing  Full  Inverse 


Magnitude  (dB)  Spectral  Magnitude  (dB) 


CG  Spectral  Estimation  —  Ideal  Correlations 
Varying  Number  of  CG  Steps  per  Angle  Increment 


Spatial  Spectrum  Estimates 


Spatial  Spectrum  Estimates 


Mm  Var 
CG  Min  Var 


2  CG  steps  per 
,05  deg  increment 
Ideal  correlations 


in 


Min  Var 
CG  Min  Var 


3  CG  steps  per 
,05  deg  increment 
Ideal  correlations 


Arrival  Angle  (Degrees) 

12  sensor  ULA:  A  =  0,05° 


Arrival  Angie  (Degrees) 

SNR1=14  dB:  SNR2=SNR3=20  dB 


CG  Spectral  Estimation  —  Ideal  Correlations 
Varying  Number  of  CG  Steps  per  Angle  Increment 


Spatial  Spectrum  Estimates 


Spatial  Spectrum  Estimates 


V 


Min  Var 
CG  Min  Var 


3  CG  steps  per 
.05  deg  increment 
Ideal  correlations 


Min  Var 
CG  Min  Var 


4  CG  steps  per 
,05  deg  increment 
Ideal  correlations 


Arrival  Angle  (Degrees) 
12  sensor  ULA:  A  =  0.05° 


Arrival  Angle  (Degrees) 

SNR1=14  dB;  SNR2=SNR3=20  dB 


CG  Spectral  Estimation 

Varying  Number  of  CG  Steps  per  Angle  Increment 


Magnitude  (dB)  Spectral  Magnitude  (dB) 


Spectral  Magnitude  (dB)  Spectral  Magnitude  (dB) 


Magnitude  (dB)  Spectral  Magnitude  (dB) 


Spatial  Spectrum  Estimates 


Spatial  Spectrum  Estimates 


MinVar 
-  CG  MV 

5  CG  steps  per 
♦05  deg  increment 
86  sensors 
J28  snapshots 
No  covariances 


50  100 

Arrival  Angle  (Degrees) 

86  sensor  ULA 


Min  Var 


CG  MV 


4  CG  steps  per 


025  dee  increment 


86  sensors 


128  snapshots 


No  covariances 


50  100  150 

Arrival  Angle  (Degrees) 

SNR1=14  dB:  SNR2=SNR3=20  dB 


CG-FFT  Spectral  Estimation:  No  Covariance  Matrix 
Number  of  CG  Steps  versus  SNR 


Spatial  Spectrum  Estimates 


Min  Var 
CG  MV 


4  CG  steps  per 
.025  deg  increment 
86  sensors 
1 28  snapshots 
No  covariances 


Spatial  Spectrum  Estimates 


Min  Var 
CG  MV 


4  CG  steps  per 
.025  deg  increment 
86  sensors 
128  snapshots 
No  covariances 


l 


50  100  150 

Arrival  Angle  (Degrees) 

NR1=14  dB;  SNR2=SNR3=20  dB 


50  100  150 

Arrival  Angle  (Degrees) 

SNR1=4  dB;  SNR2=SNR3=10  dB 


CG-FFT  Spectral  Estimation:  No  Covariance  Matrix 
SNR,  Angular  Separation,  and  Number  of  Snapshots  Decreased 


Spectral  Magnitude  (dB)  Spectral  Magnitude  (dB) 


SNRl=-6  dB;  SNR2=SNR3=0  dB 
angle  between  two  strongest 
sources  equals  10  degrees 
(diagonal  loading) 


SNRl=-6  dB;  SNR2=SNR3=0  dB 
angle  between  two  strongest 
sources  equals  5  degrees 
(diagonal  loading) 


CG  Minimum  Variance  Spectral  Estimation 
True  MV  versus  CG  MV  versus  Bartlett 


Spatial  Spectrum  Estimates  Spatial  Spectrum  Estimates 
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